## Replication code for LFSS (TESS survey) graphs presented in paper
## Ariel White (for linked fate project with Claudine Gay, Jennifer Hochschild)
## July 2015

##This code should be run AFTER running the relevant Stata code to produce the estimates plotted here.

rm(list=ls())
setwd("C://Users/arwhite/Desktop/replicationdata/tessdata") #change filepath accordingly

#load in recoded data:
library(foreign)
LFdta <- read.dta("tess2_011_hochschild_r30Mar11.dta", convert.factors=TRUE)
library(car)
lf_race1 <- recode(LFdta$lf_race, "'No'=0; 'Not Very Much' = 1; 'Some' = 2; 'A Lot'=3", as.factor.result=FALSE)
lf_racebin <- recode(LFdta$lf_race, "'No'=0; 'Not Very Much' = 1; 'Some' = 1; 'A Lot'=1", as.factor.result=FALSE)
lf_class1 <- recode(LFdta$lf_class, "'No'=0; 'Not Very Much' = 1; 'Some' = 2; 'A Lot'=3", as.factor.result=FALSE)
lf_classbin<- recode(LFdta$lf_class, "'No'=0; 'Not Very Much' = 1; 'Some' = 1; 'A Lot'=1", as.factor.result=FALSE)
lf_religion1 <- recode(LFdta$lf_religion, "'No'=0; 'Not Very Much' = 1; 'Some' = 2; 'A Lot'=3", as.factor.result=FALSE)
lf_religionbin <- recode(LFdta$lf_religion, "'No'=0; 'Not Very Much' = 1; 'Some' = 1; 'A Lot'=1", as.factor.result=FALSE)
lf_gender1 <- recode(LFdta$lf_gender, "'No'=0; 'Not Very Much' = 1; 'Some' = 2; 'A Lot'=3", as.factor.result=FALSE)
lf_genderbin <- recode(LFdta$lf_gender, "'No'=0; 'Not Very Much' = 1; 'Some' = 1; 'A Lot'=1", as.factor.result=FALSE)

LF_data <- cbind(LFdta, lf_race1, lf_racebin, lf_class1, lf_classbin, lf_religion1, lf_religionbin, lf_gender1, lf_genderbin)
head(LF_data)


#### Fig 1
library(survey)
colnames(LF_data)
mean(na.omit(LF_data$lf_race1))

##Drop group C (as discussed in paper)
nrow(LF_data)
LF_data <- subset(LF_data, LF_data$groupC==0); nrow(LF_data)
##Use weight 2 here, since I'm not comparing across experimental groups (per communication from KN/GfK).
LFdesign <- svydesign(id = ~caseid, weights = ~weight2b, data = LF_data)

##calc means for each group/var:
lfracebyeth <- svyby(~lf_race1, ~dov_ethnicity, LFdesign, svymean, na.rm=TRUE)
lfclassbyeth<- svyby(~lf_class1, ~dov_ethnicity, LFdesign, svymean, na.rm=TRUE)
lfreligionbyeth<- svyby(~lf_religion1, ~dov_ethnicity, LFdesign, svymean, na.rm=TRUE)
lfgenderbyeth<- svyby(~lf_gender1, ~dov_ethnicity, LFdesign, svymean, na.rm=TRUE)

#and also for binary (collapsed) verions:
lfracebyeth_b <- svyby(~lf_racebin, ~dov_ethnicity, LFdesign, svymean, na.rm=TRUE)
lfclassbyeth_b<- svyby(~lf_classbin, ~dov_ethnicity, LFdesign, svymean, na.rm=TRUE)
lfreligionbyeth_b<- svyby(~lf_religionbin, ~dov_ethnicity, LFdesign, svymean, na.rm=TRUE)
lfgenderbyeth_b<- svyby(~lf_genderbin, ~dov_ethnicity, LFdesign, svymean, na.rm=TRUE)


##calc means for each group/var:
lfracebyeth <- svyby(~lf_race1, ~dov_ethnicity, LFdesign, svymean, na.rm=TRUE)
lfclassbyeth<- svyby(~lf_class1, ~dov_ethnicity, LFdesign, svymean, na.rm=TRUE)
lfreligionbyeth<- svyby(~lf_religion1, ~dov_ethnicity, LFdesign, svymean, na.rm=TRUE)
lfgenderbyeth<- svyby(~lf_gender1, ~dov_ethnicity, LFdesign, svymean, na.rm=TRUE)

palette <- c(gray(.2), gray(.4), gray(.6), gray(.8))

##redo this with mean of full LF item instead of binary:
pdf("LF_figure1replication_full.pdf")
par(mfrow=c(2,2))
barplot(lfracebyeth$lf_race, main="Racial Linked Fate", names.arg=lfracebyeth$dov_ethnicity, cex.names=.55, cex.main=1.2, ylim=c(0,3), ylab="Mean reported linked fate", col=palette)
barplot(lfclassbyeth$lf_class, main="Class Linked Fate", names.arg=lfclassbyeth$dov_ethnicity, cex.names=.55, cex.main=1.2, ylim=c(0,3), ylab="Mean reported linked fate", col=palette)
barplot(lfreligionbyeth$lf_religion, main="Religion Linked Fate", names.arg=lfreligionbyeth$dov_ethnicity, cex.names=.55, cex.main=1.2, ylim=c(0,3), ylab="Mean reported linked fate", col=palette)
barplot(lfgenderbyeth$lf_gender, main="Gender Linked Fate", names.arg=lfgenderbyeth$dov_ethnicity, cex.names=.55, cex.main=1.2, ylim=c(0,3), ylab="Mean reported linked fate", col=palette)
dev.off()
# at some point I must have had a version of this code that tilted labels and added mtext, but this replicates the bars of fig1.


#################################################################################################################
## Figure 2

diffs <- read.table("sociologicaldifferences_unweighted_feb2014.csv") #load in estimates from Stata
CIlow <- diffs$diff_est- 1.96*diffs$SE  #generate 95% CIs
CIhigh <- 1.96*diffs$SE + diffs$diff_est

diffs1 <- cbind(diffs, CIlow, CIhigh)

##some plot setup:
ruler <- seq(1.5,0, by=-.5)
diffs2 <- cbind(diffs1, ruler)
diffsall <- diffs2[17:20,]
labels <- c("Upper - lower class", "Male - female", "Christian - nonChristian", "White - nonwhite")

##PLOT
pdf("sociologicaldiffs_feb2014_unweighted.pdf")
par(mai=c(.35,.15,.275,.15), oma=c(3.4,8.45,1.55,.2))
plot(NULL, NULL, main = "Linked Fate in Privileged and Non-privileged Groups", xlim = c(-1.5, 1.5), ylim = c(-.25, 1.75), xlab="Difference of means", yaxt='n', ylab="", cex.axis=.75)
axis(2, at=ruler, labels =labels, las=1, tick=FALSE, cex.axis=.95)
abline(h= .5, lty=3, col="lightgray")
abline(h= 1, lty=3, col="lightgray")
abline(h=1.5, lty=3, col="lightgray")
abline(v=0)

points(diffsall[1:4,1], ruler, pch=19) #just the first three lines: white respondents
segments(diffsall[1:4,3], ruler, diffsall[1:4,4], ruler) #CI's in the last two columns.

dev.off()



#########################################################################################################################################################
### Figure 3

rm(list=ls())

library(foreign)
coeffs1 <- read.csv("politicizationcoeffs_new.csv", header=TRUE) #load in estimates
CIlow <- coeffs1$coeff - coeffs1$serr*1.96
CIhigh <- coeffs1$coeff + coeffs1$serr*1.96 ##construct 95% CI's from the stata SE's
sig <- ifelse(CIlow>0,1,ifelse(CIhigh<0, 1, 0))#create dummy variable: does CI include 0?
coeffs <- cbind(coeffs1, CIlow, CIhigh, sig)
outcome1 <- rep(c(3, 4, 5, 6, 1, 2), 16) #reordering the outcome vars after the fact
coeffs <- cbind(coeffs, outcome1)
coeffs


pdf(file="politicizationgridplot_newlpm.pdf", width=15, height=10)

par(mfrow=c(4,6), oma=c(5,4,4,1)) #create grid of plots
op <- par(mar =c(2,1.5,1.5,.5) + 0.1) #set margins smaller to get rid of white space

for (i in 1:4){  ###this loop is for type of LF: race=1, class 2, gender 3, religion 4
for (j in 1:6){ ##this is outcome: voterreg, party, etc.

ruler <- seq(2,.5, by=-.5)
labels <- c("W", "B", "H", "A")
plot(NULL, NULL, ylim=c(0, 2.5), xlim=c(-0.750, .75), xlab="", ylab="", yaxt='n', cex.axis=1.45)
abline(v=0)
LF <- i
OUT <- j
plotco1 <- subset(coeffs, coeffs$lf==LF & coeffs$outcome1==OUT) #subset to this plot
plotco <- cbind(plotco1, ruler)
text(plotco$coeff, ruler, labels, cex=2)
##also need to make letters bold if coefficient is sig.
plotco2 <- subset(plotco, sig==1)
points(plotco2$coeff, plotco2$ruler, pch=1, cex=6)
}}
##Now label the whole plot (probably are better ways to do this)
mtext("Predicting Political Outcomes with Linked Fate", outer=TRUE,  side=1, line=2.5, font=1, cex=2.5)

mtext("Race", outer=TRUE,  side=2, line=.2, font=1, cex=1.5, adj=.9)
mtext("Class", outer=TRUE,  side=2, line=.2, font=1, cex=1.5, adj=.64)
mtext("Gender", outer=TRUE,  side=2, line=.2, font=1, cex=1.5, adj=.37)
mtext("Religion", outer=TRUE,  side=2, line=.2, font=1, cex=1.5, adj=.09)

mtext("Party ID", outer=TRUE,  side=3, line=.2, font=1, cex=1.5, adj=0.06)
mtext("Ideology", outer=TRUE,  side=3, line=.2, font=1, cex=1.5, adj=0.23)
mtext("Voter Registration", outer=TRUE,  side=3, line=.2, font=1, cex=1.5, adj=0.407)
mtext("Neighborhood Assn.", outer=TRUE,  side=3, line=.2, font=1, cex=1.5, adj=0.606)
mtext("Community Work", outer=TRUE,  side=3, line=.2, font=1, cex=1.5, adj=0.79)
mtext("Community Board", outer=TRUE,  side=3, line=.2, font=1, cex=1.5, adj=0.982)

dev.off()



